Detecting Seasonality in population indicators
#load(paste0(IO$output_kindara,"pop_indicators.Rdata"), verbose = TRUE)
load(paste0(IO$tmp_kindara, "sk1.Rdata"), verbose = TRUE)
## Loading objects:
## sk1
g = ggplot(sk1, aes(x = date, y = tot_sex)) + geom_line()
g
s = reshape(sk1, varying = list(4:ncol(sk1)), idvar = c("date","n_cycles","n_obs"), direction = "long")
s$time = colnames(sk1)[4:ncol(sk1)][s$time]
colnames(s) = c("date","n_cycles","n_obs","indicator","value")
g = ggplot(s[grep("sex",s$indicator),], aes(x = date, y = value/n_cycles, col = indicator))
g + geom_line() + ggtitle("Sexual frequency (#/number of cycles)")
g = ggplot(s[grep("sex",s$indicator),], aes(x = date, y = value/n_obs, col = indicator))
g + geom_line() + ggtitle("Sexual frequency (#/number of observations on that day)")
s$n_indicator = s$value/s$n_obs
s$n_indicator = s$value/s$n_cycles
signal = s$n_indicator[s$indicator == "tot_sex"]
time = par$time_num - 1
plot(time, signal, type = "l")
seasonal_analysis = detect_seasonal_pattern(time = time, signal = signal, title = "tot_sex Kindara", remove_weekly_patterns = FALSE)
seasonal_analysis = detect_seasonal_pattern(time = time, signal = signal, title = "tot_sex Kindara", remove_weekly_patterns = TRUE)
### TRYING WITH STL (Seasonal Decomposition of Time Series)
## weekly + yearly
# weekly period
signal_ts = ts(signal, frequency = 7)
stl_w = stl(signal_ts, s.window = "periodic")
plot(stl_w)
# yearly period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/3)
stl_y = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_y)
# bi-annual period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/6)
stl_ba = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_ba)
## "detrend" + weekly + yearly
# yearly period >> to detrend
signal_ts = ts(signal, frequency = length(signal)/3)
stl_d = stl(signal_ts, s.window = "periodic")
plot(stl_d)
# weekly period
signal_detrended = ts(apply((stl_d$time.series[,2:3]),1,sum), frequency = 7)
stl_w = stl(signal_detrended, s.window = "periodic")
plot(stl_w)
# yearly period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/3)
stl_y = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_y)
### TRYING WITH Decompose
signal_ts = ts(signal, frequency = length(signal)/3)
decompose_y = decompose(signal_ts)
plot(decompose_y)
